% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function h = plot_lambda(grid,base,lambda,N)
% h = plot_lambda(grid,base,lambda,N)
%
% Plot the hybrid variable lambda.

 switch grid.d

  case 1
   
   h = figure;

   X = zeros(grid.ns,1);
   L = zeros(grid.ns,1);
   for i=1:grid.ns
     X(i) = grid.v(grid.s(i).iv).x;
     L(i) = lambda(i);
   end
    
   [X,idx] = sort(X);
   plot(X,L(idx),'ok-');

  case 2

   csi_plot = linspace(0,1,N);
   % evaluate the basis on the csi_plot nodes
   ETA = zeros(base.nk,N);
   for i = 1:base.nk
     e_si = base.e_s{i};
     for j = 1:size(e_si,2)
       if(size(e_si,1)==1)
         ETA(i,:) = ETA(i,:) + e_si(1,j);
       else
         ETA(i,:) = ETA(i,:) + e_si(1,j)*(csi_plot.^e_si(2,j));
       endif
     end
   end
   
   h = figure;
   
   X = zeros(grid.ns,N);
   Y = zeros(grid.ns,N);
   L = zeros(grid.ns,N);
   
   for is = 1:grid.ns
      
     a = grid.v(grid.s(is).iv(1)).x(1);
     b = grid.v(grid.s(is).iv(2)).x(1);
     XY_plot(1,:) = (b-a)*csi_plot + a;
     a = grid.v(grid.s(is).iv(1)).x(2);
     b = grid.v(grid.s(is).iv(2)).x(2);
     XY_plot(2,:) = (b-a)*csi_plot + a;
      
     lambda_side = lambda((is-1)*base.nk+1 : is*base.nk);
     lambda_plot = lambda_side' * ETA;
   
     X(is,:) = XY_plot(1,:);
     Y(is,:) = XY_plot(2,:);
     L(is,:) = lambda_plot;
     
   end
   
   %plot3(X',Y',L',["k" ";;"]);
   is = 1;
   plot3(X(is,:),Y(is,:),L(is,:),'-k');
   hold on
   for is = 2:grid.ns
     plot3(X(is,:),Y(is,:),L(is,:),'-k');
   end
   
   hold off

  otherwise
   error("plot_lambda can only plot 1D and 2D solutions");
 endswitch

 title('\lambda_h');

return

